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Given a quantum Hamiltonian and its evolution time, the corresponding unitary evolution opera¬ 
tor can be constructed in many different ways, corresponding to different trajectories between the 
desired end-points. A choice among these trajectories can then be made to obtain the best com¬ 
putational complexity and control over errors. As an explicit example, Grover’s quantum search 
algorithm is described as a Hamiltonian evolution problem. It is shown that the computational 
complexity has a power-law dependence on error when a straightforward Lie-Trotter discretisa¬ 
tion formula is used, and it becomes logarithmic in error when reflection operators are used. The 
exponential change in error control is striking, and can be used to improve many importance sam¬ 
pling methods. The key concept is to make the evolution steps as large as possible while obeying 
the constraints of the problem. In particular, we can understand why overrelaxation algorithms 
are superior to small step size algorithms. 
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Classical computer simulations of quantum systems are not efficient—well-known examples 
range from the Hubbard model to lattice QCD—and Feynman argued that quantum simulations 
would do far better [|^. The essence of the argument is that quantum simulations sum multiple 
evolutionary paths (in superposition) contributing to a quantum process at one go, while classi¬ 
cal simulations evaluate these paths one by one. Formalisation of this advantage for simulating 
physical Hamiltonians, in terms of computational complexity, has improved over the years step by 
step Here, treating Grover’s quantum search algorithm as a Hamiltonian evolution 

problem, we expose the physical reasons behind the improvement in computational complexity. 

Computational complexity of a problem is a measure of the resources needed to solve it. 
Conventionally, the computational complexity of a decision problem is specified in terms of the 
size of its input, noting that the size of its output is only one bit. Problems with different output 
requirements are reduced to a sequence of decision problems, with gradually narrowing bounds on 
the output adding one bit of precision for every decision made. In such a scenario, the number 
of decision problems solved equals the number of output bits, and it is appropriate to specify the 
complexity of the original problem in terms of the size of its input as well as its output. Generalising 
the conventional classification, the computational algorithm then can be labeled efficient if the 
required resources are polynomial in terms of the size of both its input and its output. 

Popular importance sampling methods are not efficient according to our criterion, because the 
number of iterations needed in the computational effort has a negative power-law dependence on 
the precision e (i.e. Aiter as per the central limit theorem). On the other hand, finding zeroes 

of a function by bisection is efficient (i.e. Alter log e), and finding them by Newton’s method is 
super-efficient (i.e. Ater loglogg). 

1. Quantum Hamiltonian Simulation 

The Hamiltonian simulation problem is to evolve an initial quantum state |v^(0)) to a final 
quantum state \ y{T)), in presence of interactions specified by a Hamiltonian H{t): 



( 1 . 1 ) 


Alternatively, the problem can be defined as determination of the evolution operator U (T), without 
any mention of the initial and the final states. The norm of the difference between the simulated 


and the exact evolution operators specifies the simulation accuracy, say ||f/(r) — U (r)|| < e. 


We restrict ourselves here to Hamiltonians acting in finite A-dimensional Hilbert spaces. A 
general H{t) then be a dense A x A matrix, and there is no efficient way to simulate it. So we 
furthermore assume that H{t) the following features commonly present in physical problems: 

(1) The Hilbert space is a tensor product of many components, e.g. A = 2'” for a system of qubits. 

(2) The components have only local interactions irrespective of the size of the system, e.g. only 
nearest neighbour couplings. That makes H{t) sparse, with 0(A) non-zero elements. 

(3) H{t) is specified in terms of a finite number of functions, while the arguments of the functions 
can depend on the components, e.g. the interactions are translationally invariant. That allows H{t) 
to have a compact description, and consequently the resources needed to just write down H{t) do 
not influence the simulation complexity. 
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Such Hamiltonians can be mapped to graphs with bounded degree d, with vertices o compo¬ 
nents and edges •(-)• interactions. Their simulations can be easily parallelised—on classical com¬ 
puters, they allow SIMD simulations with domain decomposition. With these criteria, efficient 
Hamiltonian simulation algorithms are those that use resources polynomial in log(A^), d and log(£). 

1.1 Hamiltonian Decomposition 

Efficient simulation strategy for Hamiltonian evolution has two major ingredients. The first 
ingredient is to decompose the sparse Hamiltonian as a sum of non-commuting but block-diagonal 
Hermitian operators, i.e. H = Then each Hi can be easily and exactly exponentiated for 

any time evolution T, with exp(—retaining the same block-diagonal structure. Reducing the 
block size all the way to 2 x 2, the blocks become linear combinations of projection operators. 
Projection operators with only two distinct eigenvalues can be interpreted as binary query oracles. 

In general. Hi can be identified by an edge-colouring algorithm for graphs ||^, with distinct 
colours for overlapping edges. At most d+l colours are needed to efficiently colour any sparse 
graph. Identification of //, also provides a compressed labeling scheme that can be used to address 
individual blocks. The number of blocks is 0{m) = O(logA^), and they can be evolved simultane¬ 
ously, in parallel (classically) or in superposition (quantum mechanically). 

For example, even and odd edges of a linear chain provide a block-diagonal decomposition 
of the one-dimensional Laplacian operator, H = Ho + Hg. Its projection operator structure follows 
from H^ = 2Ho and H^ = 2Hg. The last bit of the position label identifies Ho and Hg. Eigenvalues 
of H are 4sin^(k/2) in terms of the lattice momentum k, while those of Ho and Hg are just 0 and 2. 

1.2 Evolution Optimisation 

Given that individual Hi can be exponentiated exactly and efficiently, their sum H can be 
approximately but efficiently exponentiated using the discrete Eie-Trotter formula: 

exp ( — ini') = exp ^ — i^HiT^ « ^ exp(—///,At)^ , n = T/At . (1.2) 

i i 

This approximation retains unitarity of the evolution, but may not preserve other properties such as 
the energy. The accuracy of the approximation is usually improved by decreasing At. This method 
has been used in classical parallel computer simulations of quantum evolution problems |^, 

In contrast, the second ingredient of efficient Hamiltonian simulation is to use as large At as 
possible. When the exponent is proportional to a projection operator, the largest At is the one that 
makes the exponential a reflection operator. Such an extreme strategy not only keeps the evolution 
accurate but also improves the algorithmic complexity from a power-law dependence on e to a 
logarithmic one. This is not obvious, and we demonstrate it next for the quantum search problem. 

2. Quantum Search as Hamiltonian Evolution 

The quantum search algorithm works in an A^-dimensional Hilbert space, whose basis vectors 
{|/)} are identified with the individual items. It takes the initial state l^) whose amplitudes are 
uniformly distributed over all the items, to the target state \t) where all but one amplitudes vanish. 

|VA(0)) = |.), Ivr(r)) = |t), K/|.)| = i/Viv, (/|t) = 5 ,. (2.1) 
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The simplest evolution schemes taking 1^) to |t) are governed by time-independent Hamilto¬ 
nians that depend only on 1^) and |t). The unitary evolution is then a rotation at a fixed rate in the 
two-dimensional subspace, formed by 1^) and |f), of the whole Hilbert space. In this subspace, let 



( 1/V^ \ 

Vy/(iV-l)/iVy ■ 


( 2 . 2 ) 


There are many evolution routes with U (r)|5) = |t), and we consider two particular cases in turn. 


2.1 Farhi-Gutmann’s and Grover’s Algorithms 


Grover based his algorithm on a physical intuition for the Hamiltonian [^, where the po¬ 
tential energy term |f)(t| attracts the wavefunction towards the target state and the kinetic energy 
term |5)(i'| diffuses the wavefunction over the whole Hilbert space. Both the terms are projection 
operators, and the time-independent Hamiltonian is 

tic= k)(‘5'| + |t)(t| + . (2.3) 

The corresponding evolution operator is (without the global phase) 

f/c(t) =exp(-m-at/v/i^ , h = (a/(A-1)/A,0,1/V^^ , (2.4) 


which is a rotation by angle It / s/N around the direction defined by h on fhe Bloch sphere. 

The (unnormalised) eigenvecfors of He are 1^) ± |t). They correspond fo fhe directions ±n, and 
biseef fhe initial and fhe largel slates. Thus a rolalion by angle n around h lakes |5)(5| to |t)(t| on 
the Bloch sphere, and the time required for the Hamiltonian search is T = {n/2)s/N [|T^]. 

Grover made an enlightened jump from this scenario, motivated by the Lie-Trotter formula. 
He exponentiated the projection operators in He to reflection operators; R = exp(zh/7rP) = 1 — 2P 
for any projection operator P. His optimal algorithm iterates the discrete evolution operator [jlTl], 

Ug = -(1 -2|.)(.|)(1 -2|t)(t|) = (1 - i)/ + 2/^^ia2 . (2.5) 

With Ug = exp(—///( jt), it corresponds to the Hamiltonian and the evolution step: 

^G= ;^(k)(^|-k)(l|) =/[|t)(t|,k)(^|] , Z= (^)’ 

It is an important non-trivial fact that Hg is the commutator of the two projection operators in He- 
On the Bloch sphere, each Ug step is a rotation by angle Ixs/N — \/N = 4sin 
around the direction hg = (0,1,0)^, taking the geodesic route from the initial to the final slate. 
Thai makes fhe number of steps required for Ihis discrele Hamiltonian search. 


Qt — 


cos ^{1/\/N) 


7 ^- 

4 


(2.7) 


2sin^^(l/ViV) 

Note that h and hg are orthogonal, so the evolution trajectories produced by rotations around 
them are completely different, as illustrated in Fig.|i|. It is only after a specific evolulion time, 
corresponding lo fhe solufion of fhe quanlum search problem, lhal fhe iwo irajeclories meel. 

To compare fhe rales of ihese Iwo Hamiltonian evolutions, we observe that He can be simu¬ 
lated by alternating small evolution steps governed by |5)(5| and |t)(f|, according to the Lie-Trotter 
formula. Then each evolution step governed by |t)(t| needs two binary queries [p^. On the other 
hand, Ug can be simulated using only one binary query per evolution step. 
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Figure 1: Evolution trajectories on the Bloch sphere for the quantum search problem, going from |s) to |f). 
The Hamiltonians He and Ha generate rotations around the directions h and ha respectively. 


2.2 Equivalent Evolutions 

Two Hamiltonian evolutions are truly equivalent, when their eorresponding unitary evolution 
operators are the same (upto a global phase). The interseetion of the two evolution trajeetories is 
then independent of the speeifie initial and final states. For the quantum seareh problem, we find 

Uc{T)=i{l-2\t){t\){UGf^ . (2.8) 


For a general evolufion time 0 < f < T, we have fhe relafion (similar fo Euler angle deeomposifion), 


Uc{t) = exp (/jSaa) (f/o)^' exp (/(^ + j3)a3 


(2.9) 


i.e. Uc{t) ean be generafed as Qt iferafions of fhe Grover operator Ug, preeeded and followed by 
phase rofafions. Here = 2|t)(t| — 1 is a known refleefion, and 


sin 


-1 


Qt = 


N-l 

N 


Sin 


I „ 71 L ^ 

- , p = - — - - tan tan 

2 ’^ 42 


2sin-'(l/\/iV) 


( 2 . 10 ) 


It is truly remarkable that Hg can be used to obtain the same evolution as He, even though the 
two Hamiltonians are entirely different in terms of their eigenveetors and eigenvalues! 


2.3 Discretised Hamiltonian Evolution Complexity 

Digitisation of eontinuous variables is neeessary for fault-tolerant eomputation with eontrol 
over bounded errors. But it also introduees diseretisation errors that must be kept within speeified 
bounds. The algorithmie error of the Lie-Trotter formula depends on At, whieh has to be ehosen so 
as to satisfy the total error bound e on U{t). For the simplest diseretisation seheme, 

exp ^ ^ HiAt^ = exp ( — iHiAt) ... exp ( — iHiAt^ x exp (— , (2.11) 

E^^^ = ^Y.if^i,Hj] + 0{At). (2.12) 

^ i<j 

For unitary operators X and Y, Cauehy-Sehwarz and triangle inequalities give, 

||X"-T"|| = ||(X-T)(X”^^-P...-PT"^^)|| <n\\X-Y\\ . (2.13) 
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So for the total evolution to remain within the error bound e, we need 

n\\exp{-iE^^\Atf)-I\\^n\\E^^^\\{Atf=t\\E^^'>\\{At)<e. (2.14) 


With exact exponentiation of the individual terms Hi, the computational cost to simulate a single 
time step At, , does not depend on At. The complexity of the Hamiltonian evolution is then 




(2.15) 


With superlinear scaling in t and power-law scaling in e, this small At scheme is not efficient. 

Grover’s optimal algorithm uses a discretisation formula where exp(—are reflection 
operators. The corresponding time step is large, i.e. Ate = for Eq.( |1.2[ ) applied to Eq.( |2.3[ ). The 
large time step introduces an error because one may jump across the target state during evolution 
instead of reaching it exactly. Qt is not an integer as defined in Eq.( 2.10 ), and needs to be replaced 
by its nearest integer approximation [Qt + in practice. Since each time step provides a rotation 
by angle a = 2sin^^(l/-\/iV), and one may miss the target state by at most half a rotation step, 
the error probability of Grover’s algorithm is bounded by sin^(a/2) = l/N, independent of the 
number of time steps. Since the preceding and following phase rotations in Eq.(2.9) are unitary 
operations, this error bound applies to Uc{t) as well. Thereafter, multiple runs of the algorithm and 
selection of the result by majority rule can rapidly reduce the error probability. With R runs, the 
error probability becomes less than 2 ^^^which can be made smaller than any prescribed 
error bound e. (In a drastic contrast, averaging the results of multiple runs would make the error 
probability smaller than 1 / (Ns/R) only.) The computational complexity of the evolution is thus 


With linear scaling in time and logarithmic scaling in e, this algorithm is efficient. 

To complete the analysis, we note that a digital computer with a finite register size also pro¬ 
duces truncation errors. With b-hit registers, the available precision is 5 = 2^*. With all functions 
approximated by accurate polynomials, and Euler angle decomposition reducing rotations about 
arbitrary axes to rotations about fixed axes, an individual Hi can be exponentiated to b-hit precision 
with 0{mb^) effort. The number of exponentiations of Hi needed for the Eie-Trotter formula is nl, 
which reduces to 2Qt for the Grover version. So with the choice nl5 = 0{e), i.e. b = 0(log(n/£)), 
the truncation error becomes negligible compared to the discretisation error. The cost of a single 
evolution step then scales as ^ = 0{m{\og{t/e))^), which is efficient. 


3. Extensions and Outlook 


It is straightforward to extend the preceding results to other Hamiltonians consisting of only 
two projection operators, e.g. the staggered Dirac operator for free fermions in any number of 
dimensions [13]. Eor Hamiltonians that are linear combinations of more than two projection opera¬ 
tors, e.g. three projection operators for the discretised Eaplacian on the graphene lattice, successive 
Hi can be added to the algorithm one by one in an inductive procedure. The resultant large At evo¬ 
lution is not exact, but it still has 0(1) success probability for a suitable choice of At. That keeps 
the overall scaling of the evolution efficient, C?(Zt||//||log(/t||//||/e)^) [||, 
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Our construction of the efficient Hamiltonian evolution algorithm relies on: (1) simplification 
of the Baker-Campbell-Hausdorff expansion for products of exponentials of projection operators, 
and (2) conversion of the results to a digital form allowing selection of the best one by majority 
rule. These algebraic properties are not specific to quantum computers; they can be incorporated 
in and readily benefit traditional classical simulations of quantum systems. In particular: 

(a) It is known that overrelaxation algorithms |]A], based on evolution steps that are reflections 
consistent with conservation laws, provide a much more efficient sampling of the configuration 
space (measured in terms of the autocorrelation time) than the small step size Metropolis algorithm. 
The analysis presented here provides an understanding of that observation. 

(b) Many physical problems with periodic patterns are solved using the fast Fourier transform. The 
block-diagonal decomposition provides a competitive real space method for solving them. 

(c) Projection operator decomposition of the Hamiltonian can be easily found for quantum Monte 
Carlo problems, and for molecular dynamics simulations using the Lie-Trotter formula in Euclidean 
time. The remaining task is to find a useful large step size with high acceptance probability. 

(d) Evaluations of functions other than exponentials may also simplify with the block-diagonal 
projection operator decomposition. A case of particular interest is the evaluation of the fermion 
determinant appearing in many physical problems. 

Work on such applications is in progress. 
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